Write a program to extract AVHRR sensor images from files containing recording of HRPT downlink stream from NOAA-18/19 satellites.
These files contain sequence of HRPT 10bit words stored in least significant bits of 16bit integers (using little-endian integer format). Ramaining 6 MSB are filled with 0. So the file is a sequence of 16 bit integers, one number per HRPT word.
You should decode and extract all available AVHRR image lines, select any 3 channels and save output RGB image in common image format (JPEG or PNG).
NOAA_RGB_sample.jpg
On NOAA KLM, TIP and AMSU data are updated at the major frame rate. That is, the three minor frames which make up the major frame will contain TIP data in the first minor frame, backfill in the second minor frame, and AMSU data from the AIP, in the third minor frame.
In [154]:
import numpy as np
from pathlib import Path
import struct
from matplotlib import pyplot as plt
from matplotlib.pyplot import figure
from skimage import exposure, data, filters, io, img_as_float
import cv2
from io import BytesIO
from IPython.display import Image
from ipywidgets import interact, interactive, fixed
import matplotlib as mpl
H - unsigned short (0, 255) - integer - 2 bytes
In [2]:
data = Path("./NOAA_HRPT/20180715_015444_NOAA19.raw16").read_bytes() # Python 3.5+
In [141]:
bands = 5
pixels_in_row = 2048
words_in_frame = 11090
earth_data = 10240
words_total = int(len(data)/2)
frames = int(words_total/words_in_frame)
print("frames total:", int(frames))
print("words total:", words_total)
print("check:", int(words_in_frame * frames))
print("check2:", frames * pixels_in_row * bands + (words_in_frame - earth_data)*frames)
In [6]:
arr = np.zeros(words_total, dtype=np.uint16)
for i in range(words_total):
arr[i] = struct.unpack('<H', data[2*i: 2*i + 2])[0]
In [7]:
print("First 6 symbols are: \n", arr[0:6])
print(" ")
print("They are exactly Frame Sync words:")
for i in arr[0:6]:
print("{0:b}".format(i))
In [8]:
full_scan = np.reshape(arr, (frames, words_in_frame)) # 729 frames with 11090 in each
full_scan.shape
Out[8]:
In [9]:
B1 = full_scan[0][751:10991:5]
B2 = full_scan[0][752:10991:5]
B3 = full_scan[0][753:10991:5]
B4 = full_scan[0][754:10991:5]
B5 = full_scan[0][755:10991:5]
for line in full_scan[1:]:
B1 = np.vstack((B1, line[751:10991:5]))
B2 = np.vstack((B2, line[752:10991:5]))
B3 = np.vstack((B3, line[753:10991:5]))
B4 = np.vstack((B4, line[754:10991:5]))
B5 = np.vstack((B5, line[755:10991:5]))
In [188]:
r = B1
g = B2
b = B4
rgb = np.dstack((r,g,b)) # stacks 3 h x w arrays -> h x w x 3
rgb.shape
figure(num=None, figsize=(15, 10), dpi=80, facecolor='w', edgecolor='k')
plt.imshow(rgb.astype(np.uint8))
plt.show()
mpl.image.imsave('raw.png', rgb.astype(np.uint8))
In [190]:
# Parameters
clip_limit=0.2
# Histogram adaptive equalization
rgb[:, :, 0] = exposure.equalize_adapthist(r, clip_limit=clip_limit) * 255
rgb[:, :, 1] = exposure.equalize_adapthist(g, clip_limit=clip_limit) * 255
rgb[:, :, 2] = exposure.equalize_adapthist(b, clip_limit=clip_limit) * 255
figure(num=None, figsize=(15, 10), dpi=80, facecolor='w', edgecolor='k')
plt.imshow(rgb.astype(np.uint8))
plt.show()
mpl.image.imsave('eqzd.png', rgb.astype(np.uint8))
In [166]:
# def white_balance(img):
# result = cv.cvtColor(img, cv.COLOR_BGR2LAB)
# avg_a = np.average(result[:, :, 1])
# avg_b = np.average(result[:, :, 2])
# result[:, :, 1] = result[:, :, 1] - ((avg_a - 128) * (result[:, :, 0] / 255.0) * 1.1)
# result[:, :, 2] = result[:, :, 2] - ((avg_b - 128) * (result[:, :, 0] / 255.0) * 1.1)
# result = cv.cvtColor(result, cv.COLOR_LAB2BGR)
# result = cv.cvtColor(result, cv.COLOR_BGR2RGB)
# return result
In [167]:
# figure(num=None, figsize=(15, 10), dpi=80, facecolor='w', edgecolor='k')
# img=cv2.imread('eqzd.png')
# rgb_wb = white_balance(img)
# plt.imshow(rgb_wb.astype(np.uint8))
# plt.show()
# mpl.image.imsave('rgb_wb.png', rgb.astype(np.uint8))
In [191]:
i = data.imread("eqzd.png")
In [158]:
def arr2img(arr):
"""Display a 2- or 3-d numpy array as an image."""
if arr.ndim == 2:
format, cmap = 'png', mpl.cm.gray
elif arr.ndim == 3:
format, cmap = 'jpg', None
else:
raise ValueError("Only 2- or 3-d arrays can be displayed as images.")
# Don't let matplotlib autoscale the color range so we can control overall luminosity
vmax = 255 if arr.dtype == 'uint8' else 1.0
with BytesIO() as buffer:
mpl.image.imsave(buffer, arr, format=format, cmap=cmap, vmin=0, vmax=vmax)
out = buffer.getvalue()
return Image(out)
In [192]:
arr2img(i)
Out[192]:
In [161]:
def edit_image(image, sigma=0.0, R=1.0, G=1.0, B=1.0):
new_image = filters.gaussian(image, sigma=sigma, multichannel=True)
new_image[:,:,0] = R*new_image[:,:,0]
new_image[:,:,1] = G*new_image[:,:,1]
new_image[:,:,2] = B*new_image[:,:,2]
return arr2img(new_image)
In [193]:
lims = (0.0,1.0,0.01)
interact(edit_image, image=fixed(i), sigma=(0.0,10.0,0.1), R=lims, G=lims, B=lims);